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We investigate the clustering morphology of a swarm of freely rising deformable bubbles. 
A three-dimensional Voronoi analysis enables us to quantitatively distinguish between 
two typical clustering configurations: preferential clustering and a grid-like structure. 
The bubble data is obtained from direct numerical simulations (DNS) using the front- 
tracking method. It is found that the bubble deformation, represented by the aspect 
ratio x> plays a significant role in determining which type of clustering is realized: Nearly 
spherical bubbles with x < 1.015 form a grid-like structure, while more deformed bubbles 
show preferential clustering. Remarkably, this criteria for the clustering morphology holds 
for different diameters of the bubbles, surface tension, and viscosity of the liquid in the 
studied parameter regime. The mechanism of this clustering behavior is connected to the 
amount of vorticity generated at the bubble surfaces. 



1. Introduction 

Particles dispersed in a flow can distribute inhomogeneously, showing clustering or 
preferential concentration behavior. This is attributed to the interaction between the 



two phases, and the inertia of the particles (Calzavarini et al. (2008); Toschi fc Boden 
schatz| (|2009)). A swarm of bubbles rising in a quiescent liquid is a subset of the general 
case of particles dispersed in a complex flow. This topic of bubbly flow has applications 
in bubble columns that are important in the chemical industry, in chemical processes 



such as oxidation, chlorination, in water treatment, and in the steel industry (Deen et al. 



(2000)). Freely rising bubbles in an originally still liquid are known to induce liquid ve- 
locity fluctuations which result in the so-called "pseudo-turbulence". Bubble clustering 
in pseudo-turbulence has attracted much attention because of its importance in applica- 



tions, and lack of understanding of the fundamental physics ( 


Zenit et al. (2001 


); Riboux 


et al. ( 


2010 


); Martmez-Mercado et al. (2010 


); Roghair et al. 


( 


20116)). Bunner & Tryg- 


gvason 


(2002||2003|) have conducted numerical simulations and found that deformability 



of the bubbles plays an important role in the clustering phenomenon: bubbles with small 
deformability (spherical bubbles) show a horizontal alignment, while deformed bubbles 
display a preferential clustering in the vertical direction. Meanwhile, experiments have 
found both horizontal and vertical clustering depending on parameters like bubble defor- 



mation, size, and other flow properties (Zenit et al. (2001); Cartellier & Riviere (2001); 



Martmez-Mercado et o£| ( |2010| )). 

In the present work, we revisit the issue of bubble clustering, using a Voronoi analysis 
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technique, which has been proven to be a powerful tool for quantifying the clustering 
behavior of bubbles and particles in fluid flow, see e.g. |Monchaux et al. ( 2Q1Q[ ); |Tagawa| 
et al. (2012), or Fiabane et al. (2012). Here, we extend the Voronoi analysis to study the 



geometric morphology of the clusters formed by freely rising deformable bubbles. The 
bubble data are obtained from direct numerical simulations of a swarm of rising bubbles. 

2. Voronoi analysis for clustering morphology of bubbles 

In the method of Voronoi tesselations, each Voronoi cell is defined at a particle loca- 



tion based on its neighbors (Okabe et al. (2000)). Every point inside a Voronoi cell is the 
nearest to the particle location compared to the neighbors; the exceptions being border- 
lines, vertices, and facets, which have the same distance between two or more particles. 
In a given three-dimensional distribution of particles, if the volume of the Voronoi cells 
is smaller compared to the cells in neighboring regions, then the particles belong to a 
clustering region. It has been found that a T-distribution can well describe the Probabil- 
ity Density Functions (PDF) of the Voronoi volumes of randomly distributed particles 
in three-dimensions (Ferenc & Neda (2007)), namely 



/(*) = 



3125 



l exp(— 5x), 



(2.1) 



where x = V/V is the Voronoi volume V normalized by the mean volume V. Such a 
random distribution of particles, and their corresponding Voronoi cells are shown in 
the upper panel of figure [TJb). In the lower panel of figure [TJb), the corresponding re- 
distribution fitted PDF is shown. Particles which are not randomly distributed will have 
a PDF that deviates from this T-distribution. 

Figure flta) and[TJc) show examples of preferential clustering and grid-like distribution. 
In figure Ilia), the particles prefer to aggregate in a small central region, accompanied 
by void regions. We refer to this situation as 'preferential clustering'. In this case, the 
probabilities of small and large Voronoi volumes are higher than the T-distribution. 
In figure [TJc) particles keep the same distance between each other, having the same 
size Voronoi cells. Therefore, the size distribution becomes narrower compared to the 
case of randomly distributed particles. We refer to this situation as 'grid-like structure'. 
Tagawa et al. (2012) found that these distributions can be well fitted by a T— distribution 



with a single fitting parameter a, which is the standard deviation of Voronoi volumes. 
Furthermore, this parameter a can be used for a quantification of the particle clustering. 
In this work, we use a to investigate the morphology of the bubble clustering. Here 
a is normalized by the standard deviation of randomly distributed particles (Jrnd- 

The 

indicator is (see figure [T]): cr/a rn d > 1 for preferential clustering, a/a rn d = 1 for a random 
distribution, and a /(J rn d < 1 for a grid-like structure. 

In the application of the Voronoi analysis on the present numerical data, there are two 
specific issues that are addressed below. First, the Voronoi cells of particles located near 
the edges of the domain are not well-defined, i.e. the Voronoi cells either do not close or 
close at points outside of the domain. These Voronoi cells located near the domain edges 
are usually discarded from the Voronoi analysis (as in Tagawa et al. ( 2012[ )). However, in 
the present data, the number of bubbles in the domain are small. In this case, we cannot 
afford to ignore the edge cells, as it will result in poor statistics. We take advantage of 
the periodic boundary condition of the numerics to overcome this problem. The periodic 
boundary condition enables us to form a box (3x3x3 larger, and including all particles) 
surrounding the original box. We then apply the three-dimensional Voronoi tesselation 
on the particle positions in this larger box. We can now ignore the cell edges on this larger 
box, as we have sufficient number of particles for good statistics. Also, if one considers 
the central box, although some Voronoi cells protrude into the neighboring boxes, the 
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Random distribution 
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Figure 1: Examples of the different types by which a fixed number (in this case 100) 
of particles can be spatially distributed : (a) Preferential clustering, (b) randomly dis- 
tributed, and (c) grid-like structure. The upper figures show the Vorono'i tesselations 
based on the particle positions in two dimensions, for ease of illustration. The lower 
figures show the corresponding Probability Density Functions (PDFs) of the Voronoi 
Volumes (the 3D case). The PDF corresponding to the upper figures (thick red line) are 
compared with the PDF of the randomly distributed particles (dashed black line). The 
value of the clustering indicator (i.e. the standard deviation of the PDFs normalized by 
that of randomly distributed particles, cr/a rn d) are also shown below the PDFs. 

total volume is still conserved owing to the periodic boundary conditions. This is an 
added advantage of this method. 

Secondly, the number of particles available for the Voronoi analysis is a key parameter 
that can significantly affect the results (see Fig. 6 in iTagawa et al.|fl2012| )). We check the 
dependence of the number of particles on the standard deviation of Voronoi volumes in 
figure |5|a) for randomly distributed point-like particles. We vary the number of particles 
inside each of the boxes that are replicated to form the larger box as mentioned above. 
The Voronoi tesselations are applied and the standard deviation of Voronoi volumes for 
each case of the varying number of particles is shown in figure [2j a). Also, the standard 
deviation of Vo ronoi volumes is norma lized by that of randomly distributed particles with 
numbers >10 6 (Ferenc & Neda (20071). Each error bar has been calculated by repeating 
this procedure more than 10 4 times. We see in figure (2^a) that when the particle number 
is less than 100, the value of the standard deviation changes quite significantly. We note 
the peculiarity that when a box includes just one or two particle(s), the Voronoi volumes 
are the same or half of the volume of the domain, respectively, and hence the standard 
deviation is zero. For a larger number of particles, the standard deviation grows with 
increasing number of particles and shows an asymptotic behavior, and saturates to the 
value of unity when the particle numbers approach ~ 1000. In previous work (Tagawa 
et al. (2012)) we have used a value of <r/cr r =l for the Voronoi analysis as we had 1000 
particles in our simulations. In the present work, the number of bubbles used in the 
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r-fit (eq. 2.2) 
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Figure 2: (a) The standard deviation of Vorono'i volumes as a function of the number of 
randomly distributed point-like particles in a periodic box. The standard deviations are 
normalized by the standard deviation a-p — 0-4472 for an infinite number of particles (> 
10 6 ) in a box, as shown by Ferenc & Neda (2007). In the present datasets, we consider 
16 bubbles in a periodic box, and the corresponding datapoint is indicated using a red 
triangle marker, (b) The PDFs of Vorono'i volumes for randomly distributed spheres at 
different sphere-domain length ratios D/L. The r-fit for 16 spheres expressed by equa- 



tion (2.2) with (J rn d,pp (thin black line) and the PDF for randomly distributed point-like 
particles agree well. The shape of the PDFs becomes narrower with increasing D/L due 
to the finite-size effect, (c) The standard deviations of Vorono'i volumes cr rn d normalized 
by that for the point-particle case cr rnd pp as a function of D/L. The value <J r nd/ (J rnd,pp 
decreases with increasing D/L. 

numerical simulations is 16. In this case, figure |2|a) gives the corresponding value of the 
standard deviation as (J r nd,pp/&T = 0.82. We account for this change by using a rn d yP p 
= 0.82<Jr in the equation which describes the Vorono'i volume PDF fit using the single 
parameter a ( |Tagawa et a/.|p012[ )) : 

1 



exp 



(2.2) 



This equation indeed results in a nice fit as shown in figure |2|b) , where we plot the 
Vorono'i volume PDF for 16 randomly distributed particles (blue upper triangles) and 
the curve from equation (2.2) (thin black line). 



All the above discussions were devoted to point-like particles, but in this study we 
consider bubbles with a finite-size (D = 1 — 3.5 mm). Table[l]lists the different parameters 
used in the numerics. Thus, we need to first understand the effect of finite particle size on 
the Vorono'i volume distributions. For this, we artificially generate random positions in 
three-dimensions (see § [3| for 16 perfect spheres of diameter D and change the domain 
size L to vary the sphere-domain length ratio D/L. This sphere-domain length ratio 
D/L is related to the void fraction a by the expression: D/L = (3a /Sir) 1 / 3 . For clarity 
of presentation, we have chosen to describe the clustering results using D/L instead of 
a. In figure [2^b) we show the Vorono'i volume PDFs for the 16 randomly distributed 
spheres at different D/L. The PDF of the Vorono'i volume for the randomly distributed 
point-particles and for spheres at the small value of D/L = 0.02 show quite a similar 
behavior, i.e. the finite-size effect is then negligible. The finite-size effects become more 
significant with increasing D/L, and this is seen in the shape of the PDF. The PDFs 
become narrower with increasing D/L, implying that the bubbles are distributed more 
evenly throughout the domain. At a large value of D/L, each of the spheres occupy 
relatively larger volumes in the box, which reduces the available free-space (for other 
spheres), leading to a more constrained distribution and narrower PDF shape. 
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The standard deviations of the Voronoi volume PDFs as a function of D/L are shown 
in figure [2|c). The values o~ rn d are normalized by the standard deviation obtained from 
the T-distribution fit for randomly distributed point particles, o- rn d,p P . The indicator 
o~md/o~rnd,pp decreases monotonically with D/L, starting at 1 (at D/L = 0) and reduces 
to ~l/5 for D/L = 0.36, clearly indicating the effect of finite-size. The normalization 
of the clustering indicator C = a(D / L) / a rn d(D / L) for each case used in the discussion 
below is carried out at the same bubble-domain length ratio D/L, in order to fully focus 
on dynamical effects. The total runtime for numerical simulations is about t=2 s, and 
the Voronoi volume time series reveals that the clustering is initially transient and settles 
to a quasi-steady state after t=l s. Hence, we only consider data after t=l s from the 
starting time. The clustering results are averaged over different snapshots at intervals of 
t=0.05 s. 

3. Numerical method 

Three-dimensional direct numerical simulations (DNS) have been performed to simu- 
late bubbles rising in a swarm, using periodic boundary conditions in all directions to 
mimic an 'infinite' swarm without wall effects, similarly to what has been done by |Bunner 



& Tryggvason (2002). The simulations have been carried out using a model that incor- 
porates the front-tracking (FT) method (Unverdi & Tryggvason 1 (119921)), which tracks 



the interfaces of the bubbles explicitly using Lagrangian control points distributed homo- 
geneously over the interface. Compared to interface reconstruction techniques, such as 
volume-of-fluid or level-set methods, the advantage of the FT method is that the bubbles 
are able to approach each other closely (within less than the size of 1 grid cell) and can 
even collide, while preventing (artificial) merging of the interfaces. Therefore, the size of 
the bubbles remains constant throughout the simulation. Especially for bubble swarm 
simulations with high void fractions as studied in this work, this is an important aspect. 
In addition, the interface is sharp allowing the surface tension force to act at the exact 
position of the interface. 



In our model (see Dijkhuizen et al. (2010); Roghair et al. (2011a) for details), the fluid 



flow is solved by the discretised incompressible Navier-Stokes equations on a Eulerian 
background mesh consisting of cubic computational cells: 

pd &t +pV ' ^ = ~ Vp + V,/i t W+ ( W ) T ] +^7' V -u = (3.1) 

where u is the fluid velocity and F 7 representing a singular source-term accounting for 
the surface tension force at the interface (see below). The flow field of both phases is 
resolved using a one-fluid formulation where the physical properties are determined from 
the local phase fraction; the local density p is obtained by the weighted arithmetic mean 
and the dynamic viscosity \i is obtained via the weighted harmonic mean of the kinematic 
viscosities. The interface between the gas and the liquid is tracked using Lagrangian 
control points, distributed over the interface. The control points are connected such that 
they form a mesh of triangular cells. The surface tension force is acquired by obtaining 
the pull-forces for each marker m and its neighboring cells i: F 1 ^ m = 7 (t m i x n m i). 
The shared tangent t m i is known from the control point locations, and the shared normal 
vector n m i is obtained by averaging the normals of marker m and neighboring marker i. 
Subsequently, the surface tension force is mapped to the Eulerian background grid using 



mass- weighing (Deen et al. 2004) at the position of the interface. After accounting for 
the surface tension force on all interface cells, the total pressure jump Ap of the bubble 
is obtained. The pressure jump is distributed over the bubble interface and mapped 
back to the Eulerian mesh. For interfaces with a constant curvature (i.e. spheres), the 
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Table 1: Summary of the simulation parameters. 



Case 
# 


Diameter 
D [mm] 


Void fraction 

a [%] 


Sphere- domain ratio 
D/L 


Kinematic viscosity 
v [xlO -6 m 2 /s] 


Surface tension 
7 [mN/m] 


1-3 


1.0 


10, 25, 40 


0.23, 0.31 ,0.36 


1 


73 


4 


1.0 


10 


0.23 


5 


73 


5 


1.0 


10 


0.23 


1 


7.3 


6-11 


2.0 


5 - 40 


0.18 - 0.36 


1 


73 


12-18 


2.5 


5 - 40 


0.18 - 0.36 


1 


73 


19-26 


3.0 


5 - 40 


0.18 - 0.36 


1 


73 


27-33 


3.5 


5 - 40 


0.18 - 0.36 


1 


73 



pressure jump and surface tension cancel each other out exactly on each marker, but if 
the curvature varies over the interface (which is the case for deformed bubbles), a small 
net force will be transmitted. 

At each time step, after solving the fluid flow equations, the Lagrangian control points 
are advected with the interpolated flow velocity. Spatial interpolation of the flow field 
to the control point positions is performed by a piecewise cubic spline, and temporal 
integration is performed by Runge-Kutta time stepping. Since the control points may 
move away from or towards each other, the interface mesh is remeshed afterwards, in 
order to keep the control points equally distributed on the interface (while keeping the 
volume enclosed by the dispersed element constants). 

Each bubble was tracked individually, using the locations of the control points on the 
interface to acquire the bubble position (center of mass) and bubble shape (aspect ratio), 
which were stored for further analysis. The aspect ratio is calculated from the ratio be- 
tween the major and the minor axis x along the Cartesian axes: x — (\fd x d y )/d z . Note 
that this procedure neglects diagonal shape deformations, so that strongly deformed bub- 
bles oriented diagonally may be attributed an aspect ratio of x ~ 1 (nearly spherical) . In 
this work we consider bubbles with limited deformability and under mild flow conditions 
so that these effects can safely be neglected. 

Simulations were performed using 16 bubbles in a periodic domain. For ellipsoidal 
bubbles, Bunner & Tryggvason] ( |2002| ) have indicated that 12 bubbles is the minimum 
number of bubbles that are required to simulate bubbles rising in a swarm, based on their 
terminal rise velocity. The void fraction a = Vbubbies /Vdomain was varied from dilute (a = 
0.05) up to dense void fractions of a = 0.4 by changing the domain size. In all simulations, 
the spatial resolution was determined by the bubble diameter 1.0 • 10 -3 < db < 3.5 • 10 -3 
such that the length of a cubic grid cell Ax = db/20. The physical properties represent 
typical air bubbles in water conditions, i.e. a density ratio of pliquid « 1000, dynamic 

Pgas 

viscosity ratio ^ iqmd ~ 50 and a surface tension coefficient 7= 0.073 N/m. From the 
simulations, an initial transient time of 0.2 s is discarded. The simulation parameters are 
summarized in table [T] 

An initial structured configuration of the bubble positions can cause streaming (|Bun- 



ner & Tryggvason (2003)) and liquid channeling, and may take a significant part of 



the simulation time to break up into a non-structured configuration. To prevent such 
initialization effects from influencing our simulation results, the bubble positions were 
initially set to a non-ordered fashion in the domain. Especially at higher void fractions it 
is not efficient to subsequently place a bubble randomly in the domain without allowing 
overlap. Therefore, a Monte- Carlo simulation procedure has been used to generate the 
initial positions of the bubbles which works for all void fractions (Frenkel & Smit (2002): 
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(b) Case 3 



D 

DIL 



1.0 mm 
0.23 
y =73 mN/m 
C =0.62<1 

(Grid-like structure) 




D ~- 
DIL -- 



1.0 mm 
0.36 
y =73 mN/m 
C = 0.36 < 1 

(Grid-like structure) 



D ~- 
DIL -- 



1.0 mm 
0.23 

) y =7.3 mN/m 
C = 1.22 > 1 

(Preferential clustering) 



(d) Case 33 




3.5 mm 
0.36 
y =73 mN/m 
C =1.12>1 

(Preferential clustering) 



Figure 3: Snapshots of bubbles showing typical clustering morphologies. The periodic box 
is indicated by the thin black lines. The bubbles are colored to aid readers to distinguish 
between individual bubbles. 



Beetstra ( 2QQ5[ )). First, the bubbles are placed as spheres in a structured configuration in 
the domain in a simple cubic configuration. Depending on the void fraction, the bubbles 
might overlap with each other. We now define the potential energy of the system as: 
E = [\xi — Xj\/(Ri + Rj)] n with a variable n characterizing how steep the potential is. 
The position of a bubble i is given by Si and its radius by Ri. Each bubble is now moved 
by a small amount in a random direction and the potential energy is determined. A move 
is accepted whenever the potential energy remains equal or becomes less, whereas a move 
that increases the potential energy is accepted only if it is smaller than a critical number 
c: c = exp[/c(^oid — ^new)], where k was set to 50. In a single iteration, each bubble is 
allowed 200 (attempts of) displacement. Then, the power n is gradually increased from 
an initial value of 6 up to a final value of 100, and a new iteration starts. The potential 
energy of the system as a whole decreases during this process, and when the final state 
has been reached and the bubbles show no overlap at all, the positions of the bubbles are 
accepted for use as starting positions in the front tracking model. Additionally, in the ini- 
tial transient of the front tracking simulations, the bubbles accelerate, deform and move 
through the domain, which also changes their relative positions. This start-up stage is 
discarded from further analysis. For the random positions of the point-particles as shown 
in figure [2jc), the same procedure was used, except that we chose the number of allowed 
displacements for each particle per iteration to be 10 4 . 

4. Results and Discussions 

We recall that the value of the clustering indicator C = cr/a rn d are C > 1 for preferential 
clustering, C = 1 for a random distribution, and C < 1 for a grid-like structure (figure [lj. 
In figure [3j we now show the typical bubble clustering snapshots. Figure |3|a) displays 
a snapshot of the case of bubble diameter D = 1.0 mm at D/L = 0.23. We observe 
horizontal clustering in one layer. This snapshot corresponds to C < 1, a grid-like cluster 
morphology. In figure [3jb), D = 1.0 mm at D/L = 0.36, and the bubbles show horizontal 
clustering in a double layer, owing to larger D/L (i.e. larger void fraction a), and the value 
of C is correspondingly less than 1. It must be noted that the shapes of the bubbles for 
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Figure 4: Left panel: (a) Bubble clustering results: A contour plot of the clustering in- 
dicator C as function of the bubble-domain length ratio D/L for all bubble sizes D (in 
mm). The colorbar indicates the magnitude of C. All the simulation cases in table [I] 
are shown here except 4 and 5. Right panel: The normalized standard deviation of the 
Vorono'i volumes C, as a function of the aspect ratio x f° r: (b) fixed bubble size, D = 
1.0 mm (cases 1, 2, 4, 5 in table |T|, and (c) fixed bubble-domain length ratio, D/L = 
0.36 (cases 3, 18, 26, 33 in table]!]). The standard deviation increases with increasing 
aspect ratio, indicating that the shape of the bubbles plays a crucial role in determining 
the clustering morphology. Spherical bubbles with an aspect ratio \ ^ 1-015 have C < 1, 
indicating the grid-like structure. All deformed bubbles with x ^ 1.015 have C > 1, 
implying preferential clustering. 



D = 1.0 mm at (a) D/L = 0.23 and (b) D/L = 0.36 are almost spherical. In figure [3]^c), 
we show a case of D = 1.0 mm at D/L = 0.23 with lower surface tension, where the 
bubbles are evenly distributed and the horizontal one-layer clustering no longer prevails. 
In figure |3|d), in the case of D = 3.5 mm at D/L = 0.36, the bubbles are more evenly 
distributed. The corresponding C values for both cases are larger than 1, i.e. indicating 
preferential clustering. The bubbles with D = 1.0 mm at D/L = 0.23 and those with D 
= 3.5 mm at D/L = 0.36 have a deformed shape. 

For a quantitative discussion, in figure [4^ a) we show (as a color contour plot) the values 
of the clustering indicator C at different bubble-domain length ratios D/L and bubble 
sizes D. The data points in the parameter space (table [T]) are indicated using open blue 
circles. The formation of grid-like clusters (C <1) were only encountered in the cases of 
1.0 mm diameter bubbles at D/L = 0.23 and 0.36 (figure [4^ a)). All other cases show 
preferential clustering (C >1), but to a different extent. 

It is well-known that a rising spherical bubble with a free-slip boundary condition 
generates little vorticity, whereas a rising deformed bubble has a wide region of wake 
structure behind it (Magnaudet & Mougin (2007)). The amount of vorticity generated 
from the bubbles determines the clustering morphology. The flow around spherical bub- 
bles can be expected to be close to potential flow, containing little vorticity, and these 
bubbles form a grid-like cluster (in the horizontal plane). Deformable bubbles with large 
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Figure 5: The maximum number of bubbles in a horizontal plane. The dotted line shows 
the values for square packing. The 1 mm bubbles (cases 1-5 in table [I]) show a increasing 
trend with decreasing D/L indicating that they form horizontal clusters. The 2.5mm 
bubbles (cases 12-18 in table [T]) show almost constant values of 6, implying that there 
is no horizontal clustering. The result in the case of liquid with high viscosity (case 4 in 
table [I] green triangle) and that for low surface tension (case 5 in table [I] red square) 
are also shown. 

wake regions show aggregation in the vertical direction. Hence, in the discussion of the 
results, we focus on the bubble shape (or deformability) characterized by the bubble 
aspect ratio \. Below, we fix the size and bubble-domain length ratio, and discuss the 
clustering at different x- 

First, we keep the bubble size constant (D = 1.0 mm) and discuss results for different 
bubble-domain length ratios, surface tension, and viscosity. Figure [4]^b) shows the values 
of the clustering indicator C as a function of the bubble aspect ratio x- The value of C 
increases with an increase in the aspect ratio, indicating that the bubble shape is crucial 
for the clustering structure. We now fix the bubble-domain length ratio (D/L = 0.36) 
and study the clustering behavior for different bubble sizes in figure [4^c) . Although the 
fixed parameter is different in this case, we still see the same trend of increasing C with 
increasing aspect ratio x- Overall, we find that the shape of the bubbles is crucial for 
the structure of the clustering. Spherical bubbles with x ^ 1.015 have C < 1, indicating 
the grid-like clustering structure. All the deformed bubbles with x ^ 1.015, have C > 1, 
indicating the preferential clustering morphology. 

We also quantify the intensity of horizontal clustering by counting the maximum num- 
ber of bubbles at a horizontal plane, for the cases of D = 1 mm and 2.5 mm bubbles. 
The bubbles are sliced at the horizontal plane and divided by the area of a circle based 
on the bubble radius. Figure [5] shows the maximum number of bubbles in a horizontal 
plane versus the bubble-domain length ratio. The line obtained from the theory of square 
packing is also shown for the sake of comparison. On one hand, the 1 mm bubbles show 
a trend similar to the theoretical line, indicating that they organize themselves to form 
horizontal clusters. On the other hand, the 2.5 mm bubbles show almost constant values 
around 6, indicating the absence of horizontal clustering. The result of the lower surface 
tension case (case 5 in table [I]) is rather close to the 2.5 mm case (case 13 in table [I]) 
due to deformation. These results for the horizontal clustering are consistent with those 
obtained from the Vorono'i analysis. 

5. Conclusion 

In this work, we have applied the three-dimensional Vorono'i analysis on DNS data of 
freely rising deformable bubbles in order to investigate the clustering morphology. The 
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numerics used a front-tracking method which allows the simulation of fully deformable 
interfaces of the bubbles at different diameters, bubble-domain length ratios, surface ten- 
sion, and liquid viscosity. The present Vorono'i analysis takes into account the number of 
bubbles and finite-size effects. It then provides a clustering indicator C = cr/cr rn ^, where 
a is the standard deviation of Vorono'i volumes of the bubbles and a rn d is the standard 
deviation of Vorono'i volumes of randomly distributed particles. We quantitatively iden- 
tify two different clustering morphologies: C > 1 for preferential clustering and C < 1 for 
a grid-like structure. Our results indicate that the bubble deformability, represented by 
its aspect ratio x> plays the most crucial role in determining the clustering morphology. 
The grid-like morphology is observed in the case of nearly spherical bubbles with x ^ 
1.015. When the bubbles are deformable, for x > 1.015, a preferential clustering behavior 
is observed. This clustering behavior is believed to be related to the amount of vorticity 
generated by the bubbles. The preferential clustering for deformed bubbles is due to the 
low pressure regions in their wakes, which attract other bubbles. Spherical bubbles tend 
to form a grid-like structure due to reduced vorticity generation. 
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